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ABSTRACT 


An iterative optical vector-matrix multiplier with a microprocessor-controlled 
feedback loop is capable of performing a wealth of diverse operations. In this paper, 
we survey and describe many of these operations to demonstrate the versatility and 
flexibility of this class of optical processor and its use in diverse applications. 


1 . INTRODUCTION 


The optical vector-matrix multiplier [1] is a general purpose optical processor. 
The addition of a microprocessor-controlled feedback loop results in an even more 
general purpose and far more powerful optical processor [2-6]. In this paper, we sur- 
vey and describe a selected set of the operations achievable on such a processor. 

In Section 2, a general description of the system is advanced. This is followed 
in Section 3 by a description of how bipolar and complex-valued data are handled on 
the system and how the convergence of the iterative algorithm is insured. We then 
address in Section 4 its use in the solution of linear difference and differential 
equations and linear algebraic equations. In Section 5, we consider application of 
this processor for the solution of the least-squares problem. In Section 6, we 
address its use for deconvolution and for the computation of the eigenvalues and 
eigenvectors of a matrix. In Section 7, our attention turns to the use of the system 
for matrix-matrix multiplication, the solution of linear matrix -matrix equations and 
matrix inversion. We then consider in Section 8 its use in the solution of nonlinear 
matrix equations with specific application to the linear-quadratic-regulator problem 
and algebraic Riccati equation of optimal control engineering. 


2. ITERATIVE OPTICAL PROCESSOR SYSTEM 


The general schematic of the J^terative optical p^rocessor (lOP) is shown in Fig- 
ure 1. This figure illustrates the use of the lOP both as a vector-matrix multiplier 
and as an iterative processor for the solution of linear vector-matrix equations. 
Since an understanding of this architecture is paramount to the remainder of this 
paper, we review its operations. To multiply a vector by a matrix, the elements of 
the vector are used to spatially modulate a linear array of flight emitting ^iodes 
(LEDs) or j^aser _diodes (LDs) at P^. We denote the light distribution leaving P]^ by 
the vector x. Plane P^ is imaged vertically onto P 2 and the output from each LED is 
expanded to uniformly illuminate the corresponding row of a 2-D mask at P 2 . The light 
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distribution leaving each column of P 2 is then summed onto a different photodetector 
at P 3 . With the transmittance of P 2 specified by the matrix the vector output from 
the linear photodetector array at P 3 is the matrix-vector product as described in 

[1] and earlier by others [7,8], 

In our system [2,3], we have included feedback of the photodetector outputs to 
the LED inputs through an electronic feedback system. In our original description of 
this system [2,3], the mask was [_! - H] 5 where 1 is the identity matrix. The P 3 out- 
put is then [1 - where x(j) denotes the vector x at the j-th iteration. An 

external vector was added to this matrix-vector product in the electronic feedback 
system to form the new iterative input x( j + 1) . Our system thus implemented the 
iterative algorithm 


x(j + 1 ) = U “ H]x(j) + y = x(j) - Hx (j) + Z* ( 1 ) 

In the steady-state when 5 £(j + 1) - x(j), equation (1) reduces to ^X“ Z out- 

put X is the solution 

X = (2) 

of the vector-matrix equation 

Hx = y. (3) 

In the newest [5,6] version of this lOP, we have: (1) used fiber optics to 

realize the required projection from P]^ onto P 2 (this greatly improves the system’s 
accuracy as well as its mechanical and positional stability and reduces its size and 
weight); (2) placed the photodetector at P 3 in direct contact with the matrix mask at 
P 2 (this is possible when the 4 mm height of each photodetector is matched to the 
height of the matrix mask at P 2 . This further reduces the size and weight of the sys- 
tem and makes it even more stable for airborne applications); (3) included a micro- 
processor system with an arithmetic ^ogic unit (ALU) , memory and hardwired multiplier 
in the electronic feedback loop (Figure 2), (this increases the system’s flexibility 
and versatility); and (4) modified the feedback system and the iterative algorithm to 
incorporate an acceleration parameter to insure convergence of the iterative algorithm 
(Section 3) . 

When funding permits, we plan: (1) to increase the size of the P^ input and the 

P 2 mask from its present 10 x 10 level; ( 2 ) to incorporate a real-time and reuseable 
spatial-light modulator electro-optical mask element (such as a CCD-addressed liquid 
crystal light valve [9]) into the system; and (3) to increase the repertoire of opera- 
tions and applications for the system. In Sections 4-8, we advance the first descrip- 
tion of a selected number of general operations and problems as well as applications 
for which this new optical processor can be used. 


3. CONVERGENCE AND HANDLING BIPOLAR AND COMPLEX-VALUED DATA 


In this section, we first detail how this non-coherent optical system can be 
used to operate on bipolar and complex-valued data [5,6]. Since the LED outputs at 
P^, the transmittance of the mask elements at P 2 and the photodetector outputs at P 3 
are all real and positive, and since the dynamic range of the mask at P 2 is finite, 
various scaling, biasing and data encoding techniques are required to process vectors 
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and matrices with bipolar and complex-valued data on the system. To handle complex- 
valued data, we partition the matrix H into its real and imaginary parts and H-, 

respectively, as 


H = 


H -H. 

— r — 1 


H, 

~x 


H 

—r 


(4) 


where H and H. are bipolar. This requires a P 2 mask with four times the space-band- 
width product of H, an input LED array with twice the number of elements in >£ and a 
linear detector array with twice the number of elements present in y^. 

To handle bipolar data, we decompose the input vector into its positive and 
negative components. Let _a"^ and ^ denote the positive and negative components of 
the input vector ^ used in the actual system. The M elements and a 2 ^ of ^ and 

are generated according to 


(5) 


where the are the elements of the bipolar physical vector x. 

The system is then operated twice, first with as the input and then with ^2 
as the input (with the same physical mask ^ used for all cycles) . The outputs from 
the system on successive cycles are ^_a^ and ^^ 2 * These outputs are: (1) subtracted; 
(2) scaled by (h “ h) , where h and li are the maximum and minimum values of the ele- 
ments of the matrix H ; and (3) biased ^ ^ system’s output ^ after two 

successive cycles is thus '-m=l m 




To insure that the transmittance of the physical mask jB has a transmittance for each 
element satisfying 0 ^ b^^ 1, we scale and bias ^ such that 



where h denotes the values of the elements of H and denotes the transmittances 
of the pnysical mask ]B placed at P 2 . 

The second issue to be detailed in this section is how rapid convergence of the 
iterative algorithm is insured. This is achieved by modifying our original iterative 
algorithm to include an acceleration parameter O) (as shown in Figure 1) . The new 
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system thus implements the iterative algorithm [4-6] 


x( j + 1) = x( j ) + “tz - H X ( j ) ] . 


( 8 ) 


We select 0) by one of the following criteria to insure rapid convergence of the algo- 
rithm. We can select 


0) = -1/X 


max’ 


(9) 


where A is the maximum eigenvalue, in absolute value, of the matrix H, or we can 
max ^ ^ ^ 

obtain an approximation to (.9) by [10] 


A < 
max — 


|h1 


M N 

E E h 

mn 


1/2 


( 10 ) 


where | |h( | is the Euclidean norm of the (N x N) matrix H. The O) criterion resulting 
from the upper-bound in (10) can be modified by selecting O) = “'^^/^max^ where K > 1 is 
a constant selected empirically from analysis of a specific problem. (In several spe- 
cific case studies that we have considered, K - 2-3 was found to be adequate.) The 
microprocessor feedback system (Figure 2) performs the necessary scaling, biasing and 
preprocessing described by (5) and (7), the detector post-processing in (6), and the 
acceleration parameter selection noted in (10), In Section 6, we describe how the 
lOP itself can be used to calculate the acceleration parameter 03 in (9) . 


4. SOLUTION OF SIMULTANEOUS LINEAR EQUATIONS 


As our first general lOP application, we consider the use of the system of Fig- 
ure 1 for the solution of simultaneous linear (difference, differential or algebraic) 
equations. The general iterative algorithm for the solution of linear difference 
equations is 


x(j + 1) = £x (j) + 1 > (11) 

where j is the discrete- time index, f_ is the forcing function equal to a constant vec- 
tor (or, more generally, a vector of time-functions) , x(j) is the state of the system, 
and ^ is the open-loop system matrix. 

The first application of the TOP of Figure 1 to the problem in (11) is as a dy- 
namic system simulator. In this case, we model the physical system by the state-space 
differential equation model 


dx 

— = A X + b , (12) 

dt 


where A is the open-loop system matrix and the vector ]b is a constant or function of 
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time. We then utilize a numerical analysis algorithm to discretize this system model. 
We illustrate the approach by the f orward-Euler approximation [11] 


dx(t) 


dt 


t=jT 


x(j + 1) - x(j) 
T 


(13) 


where T is the constant time- increment (or step-size) between discrete samples de- 
scribed by the index j. Substituting (13) into (12), we can simulate the physical 
system modeled by (12) on the lOP of Figure 1 by the linear iterative algorithm 


x(j + 1) = [I f TA]x(j) + Tb . (14) 

By analogy with (11), [jL ^ A A “ A • Other numerical analysis algorithms, 

such as the trapezoidal rule, are possible. The most preferable ones appear to be 
the Runge-Kutta and predictor-corrector algorithms [12]. 

The second application of (11) on the lOP is the iterative solution of linear 
algebraic equations. In this case, the iterative system algorithm of (8) is used. 

When rearranged in the form in (1) , we can identify the ^ matrix and the ^ vector in 
(11) as [I_ - 0)H] - ^ and ocy = f_. The algorithms in (8) and (1) have been used [2,13] 
to calculate the set of adaptive weights necessary in adaptive phased-array radar sig- 
nal processing. While space does not allow us to elaborate on this application, we 
note that, in this case, the matrix H in (3) is the covariance matrix M of the far field 
input to the antenna, y is the steering vector ^ for the antenna and the unknown vec- 
tor 2 ^ to be determined is the set of adaptive weights w to be applied to the received 
signals to steer the antenna in the desired direction (defined by ^) and to null the 
noise field (defined in frequency, velocity, angle and range by M) . 

Four iterative algorithms to solve (3) for k given by (2) can directly be iden- 

tified. The first is the Richardson algorithm [14] of (1) implemented with the accel- 
eration parameter o) as in (8) . The remaining three algorithms can be described in a 
new and quite useful formulation by decomposing the matrix II into the sum of a diago- 
nal matrix ^ and lower and upper triangular matrices and U as 

H = D - L - U . (15) 


In terms of the h and U, we can write the remaining three algorithms as [15]: 


Jacobi: 
Gauss-Seidel : 
Over relaxation: x(j + 


x(J + 1) = 
x(j + 1) = 

1) = j(D - 


^(L + x(j) + D V 

[(D - + (D - L)”^: 


(OL) 1[(1 - 03)D + 0)U] 


X(j) + w(D - 0)L) 


(16a) 

(16b) 


Z.- 

(16c) 


The choice of one of the four algorithms in (8) or (16) depends on many factors that 
are highly application and problem dependent (e.g., convergence of the algorithm, dy- 
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dynamic range of the mask, number of iterations required and the need for an adaptive 
mask) . 


5 > LEAST-SQUARES PROBLEMS 


The solution to many physical and statistical problems can be formulated in the 
general context of a least mean-squares problem. A model is postulated to approximate 
measured data and the parameters of the model are selected to minimize the mean-square 
error between the measured data and the model. Curve-fitting and linear and poly- 
nomial regression are classical examples of least mean-squares problems. We formulate 
and address below the solution on our lOP of this class of important problems. 

We begin by reformulating our simultaneous linear algebraic equation problem 
solution in Section 4 as a least-squares problem. in the context of Figure 3, the 
input to the operator H is x and noise or an error source ^ is present at the output such 
that the observable output y_ differs from Ux . The problem is thus to find the x 
that minimizes the square of the difference 

J = Ily - I ^ 


between y_ and H x . 


Noise or 
Error Source e 


V ^ 

H 

I 

X IP 

— i 

Hx 



FIGURE 3 Schematic diagram of 

a least-i 




Upon expanding (17), we obtain 


J = Iz - [y “ = zz ~ y'^ii2£ - ^^z + 2i^li'^Hx 


(18) 


To minimize (18) , we set the partial derivative or gradient of the scalar per- 
formance index in (17) with respect to the unknown vector x equal to zero. The re- 
sultant column vector is then the solution of 
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or 


(19) 



In least-squares data processing, y is an M-vector, H is an (M x N) matrix and x is 
an N-vector, where M ^ N, We consider first the case when the number of unknowns N 
equals the number of equations M. In this case, H is a square matrix and the solution 
to (19) becomes 


X = 


T 

H H 


T -1 -T T -1 

H y = H H H y = H y 


( 20 ) 


where H is assumed by implication to be a non-singular matrix and therefore to be in- 
vertible. This occurs when the y = Hx problem being solved is well-formulated as 
occurs in all practical engineering problems. Thus, in this case, the solution x = 
minimizes the quadratic performance index in (17). For this case, we solve the 
least mean-square problem by our Richardson algorithm of (8) . 

A more interesting problem (corresponding to the practical case of curve-fitting 
occurs when there are more equations M than unknowns N (i.e., when M > N) . In this 
case, H is non-square and therefore non- invertible. We thus solve (19) by the new 
iterative algorithm 


I 

1 

Pm! 


r p 

1 

x(j + 1) = x(j) + 0)< 

1 

1 

[h y 

- 

T 

H H 

x(j) 

i 

1 

L _ 


— 

J 


T 

This is equivalent to pre-multip lying (3) by H and applying our iterative Richardson 
algorithm in (8) to the solution of [H*^H]x = ^y rather than H x = 

The least-squares fitting of experimental observations leads to a more general 
iterative vector-matrix solution. Suppose that we have L data points p^ = {p^} and 
wish to approximate the set of observations z(p^), in the least-squares sense, by the 
finite linear combination 


zip) 




( 22 ) 


of basis functions ^ = ^^n^ with the associated weighting coefficients _c = In 

(22), each of the N basis functions ^ are evaluated at the L data points and the 

relationship in (22) is approximate if there are more data points L than basis func- 
tions N. The errors or residuals [z(p^) - ^"^^(p.) ] of the curve-fit are the differences 
between the observed data z(p^) and the approximation c^^(p^) • According to the princi- 
ple of least-squares , we select the weighting coefficients in (22) to minimize the 
sum-of-squares of the residuals . To find the coefficients c_, therefore, we minimize 
the mean-square difference 


J = 


DataL 


-(pj) - 




(23) 
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We proceed as before. Upon setting 9j/3_c = _0, we find 


E 

Data 

T 

i(P£>i (P£ 

c = 



— 


'p.J 


H 

X = 



(24) 


which is again of the general form Hx = y_^ with the summation over the {p^} data 
points incorporated into the matrix and vectors in the algorithm in (21) . When the 
basis functions ^ = ^^n^ specified, the matrix H can be precomputed and the prob- 

lem solved by the iterative algorithm in (8) or (21) . A simple discrete example of 
{<(>n> arises in the calculation of the best straight-line fit through a set of data 
points. In this case, = 1 and 4^2 ~ U ^(p) ~ ^1 ^2P* ^ simple continuous 

example case is the Fourier series expansion of z(p) , in which case the are the 

complex exponentials exp (+jn(jOQp) . 

In this section we have illustrated and formulated the least mean-squares prob- 
lem as the iterative solution of the system of linear algebraic equations Hx ~ y_ on 
the lOP. We recognize that the range of applications of our formulation includes 
curve-fitting, linear and nonlinear regression, state and parameter estimation, orbit 
determiniation and signal processing in control and communication engineering. We 
will address these applications in our future work. 


6. DECONVOLUTION AND EIGENVALUE PROBLEMS 


The need to implement a deconvolution frequently arises. In such applica- 
tions (Figure 4) , the measured output data y^ = ^^m^ from a system characterized by 
the impulse response {hj^} will be a modified version of the original input vector x = 
{x^}; i.e.. 


m 


m 

E 

n=0 


(m-n) 


X . 

n 


(25) 


In (25), we assume that the system is causal, otherwise we appropriately shift the im- 
pulse response so that “ 0 for m < 0. The solution of (25) for the unknown input 
2£ is directly described by the vector-matrix equation y = , where (for the case of 

]q = ]Sf = 3) the matrix H is shown below: 



(26) 
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= { x„ 


■y 




FIGURE 4 Simplified schematic diagram of a deconvolution processing problem. 


In the general case of M > N = 3, the size of H grows, but H still has only 
three non-zero diagonals (with equal values hg, hj^ and h 2 , respect ively^ along each 
diagonal), and there are an additional (M - N) elements with zero value added to 
The solution of the deconvolution problem in (25) is thus directly realizable on the 
lOP by the iterative algorithm in (8) with the matrix mask in the form shown in (26). 
Once again, we see how different problems can be solved on the same lOP by different 
choices of the matrix mask and the iterative algorithm used. 

Another quite general and useful matrix operation is the calculation of the 
eigenvalues and corresponding eigenvectors of the (N x N) matrix H, This opera- 
ation can be performed by the iterative algorithm [16] 


x(j + 1) = H X (j) 


(27) 


which is equivalent to (8) with the exogenous vector ^ set equal to zero. To see how 
(27) allows calculation of the and we assume that | X^ | is the largest eigen- 
value and that the eigenvalues are ordered, in absolute value, according to [X^j > 

1^2 I > •••• > |X]^|. By singular value decomposition, we write H and the original 
initialization vector 21 ^( 0 ) as 


and 


N T 

H = E 4 . X . ^7 

- i=l 


N 




(28) 


After j iterations, 2i(j + 1) = obtained. Substituting (28) into this ex- 

pression, we find (for large j) 


x(j + 1) = 


= u3 


H-'Xj^(O) 


Large j 




(29) 


We form the component-wise ratio x^(j + l)/xj^<j) and from it find X]_ (the largest 
eigenvalue of H) . We divide the denominator of the output vector x by the sum of the 
squares of its elements and thus obtain the principal eigenvector We then use the 

new initial vector ~ Repeating the same iterative procedure in 
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(27), we then find X 2 and ^ 2 * continuing this process, all of the eigenvalues and 

eigenvectors of H can be determined in decreasing order. Modif ications to this pro- 
cedure allow us to find the eigenvector whose eigenvalue lies closest to a prescribed 
value [16] . Many other extensions of these techniques and applications of the above 
results are possible and will be the subject of future work. 


7 . SOLUTION OF MATRIX EQUATIONS 


Many problems involve matrix-matrix multiplication. There are two ways to real- 
ize a matrix-matrix multiplication on our vector-matrix multiplier. In the case of 
large matrices, the preferable technique is to form the product Y = M • ^ of the 
matrices M and N by feeding sequentially the columns n^ of N as successive input vec- 
tors. We thus realize the matrix-matrix product by performing N vector-matrix multi- 
plications; i.e.. 


= M • N = 


Bn ^ 


(30) 


where the matrix output Y appears sequentially as N column vectors; i.e., Y = [yqX2 ^ 
. . y^ ] . A second technique is to vectorize the matrix N into an element vector 
whose first N elements are the elements of Un and the next N elements are the elements 

— jL o o 


of n 


2 ^' 


9 2 

The matrix M is formatted as an (N x N ) matrix whose diagonal blocks 


are replications of the matrix M; i.e,. 


M • N 



n^ 




Sn 


^1 




(31) 


With a matrix-matrix multiplier realized by (30) or (31) , we can thus use the 
lOP to solve matrix-matrix equations such as 


HX = Y 


(32) 


by vectorizing the matrix Y or by operating the system on sequential cycles with the 
column vectors of Y as successive inputs. A particularly useful operation that is 
now possible is matrix inversion. In this case, we solve (32) using (8), but with 
Y = ^ (the identity matrix) . A final useful matrix equation that frequently arises 
is the solution for the embedded matrix X in an equation of the form 


AXB = Y 


(33) 


In this case, we write the matrices IT and X as the column vectors _C[^] and _C[X], whose 
elements are the lexigraphically-ordered elements of the matrices Y and X, respective- 
ly. We solve (33) for X by writing (33) in the form 

CA © 1^]C[X] = C[Y] , (34) 
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where 

A, 



denotes the outer or Kronecker product. 


For the case of the (2 x 2) matrix 


A 






(35) 


T 

Solution of the Lyapunov matrix equation A ^ ~ X possible on the vector- 

matrix lOP. 


8 . SOLUTION OF NONLINEAR MATRIX EQUATIONS 


In Section 7, we described how we have broadened the repertoire of operations 
achievable on the lOP to include the full-class of linear and embedded matrix equa- 
tions. As our final general iterative vector matrix operation, we consider its use 
in the solution of nonlinear matrix equations. The need for the solution of such 
equations arose in our original use of the system for adaptive phased-array radar 
processing [2] and in our present optimal control linear-quadratic-regulator applica-- 
tion [4] . 

The general problem of two quadratic equations in the two unknowns p and q can 
be written as 


2 2 

f (p.q) =ap +bpq+cq +dp+eq+r =0 
n n n n n n n 


(36) 


for n = 1 and 2. We rewrite this pair of two nonlinear equations in vector form as 

f [x] - 0 (37) 

T T 

where ^[x] = [fi(x) ^2^— ” CP then solve (37) by the Newton-Raphson 

iterative algorithm 

-1 

f[x(j)]. (38) 

x(j) 


x(j + 1) = x(j) - 


Sf 


3x 


The solution of (38) requires two iterative loops. The Jacobian matrix _J = [3f_/82£] 
is stored algebraically. At each iteration, evaluated numerically with dedicated 

electronics at the last iterate * The inverse matrix [_J [2l(j) 1 evaluated on 

the optical vector-matrix system in an inner iterative loop. The new 2l(j + 1) iterate 
is then evaluated optically in an outer iterative loop. 

Three immediate applications for such a nested two-loop iterative algorithm 
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arise. The first occurs in the implementation of the overrelaxation linear algebraic 
equation solution of (16c) with - (j)L] on the left-hand side of the equation. In 
this case, a matrix-vector multiplication is required to obtain the external vector 
to be added. The second case arises in the use of (8) to solve for the adaptive 
weights w for a phased-array radar whose noise field is described by the covariance 
matrix M and the direction of its main lobe is defined by the steering vector _s. In 
this case, we solve the vector-matrix equation M • w = ^. When the noise distribution 
varies, the matrix mask M must be updated. 

The lOP application we are presently considering under support from the NASA- 
Lewis Research Center is the use of the lOP in the solution of the Mnear-jquadratic- 
r^egulator (LQR) problem of optimal control. In this application, the lOP is used to 
calculate the optimal controls to be applied to an PlOO aircraft engine. To determine 
these control signals, we must solve the algebraic Riccati equation and calculate the 
LQR feedback gains. In this application, we use: (1) the vectorization of a matrix; 

(2) the Kronecker product technique; and (3) the nested inner and outer loop system 
to solve the nonlinear algebraic Riccati equation [4]. 


9 . SUMMARY 


In this paper, the general-purpose nature of an iterative vector-matrix pro- 
cessor has been emphasized. This system is capable of solving a wealth of general 
purpose applications. General operations described included: linear difference and 

differential equations, linear algebraic equations, matrix equations, matrix inver- 
sion, nonlinear matrix equations, deconvolution and eigenvalue and eigenvector compu- 
tations, Engineering applications being addressed for these different operations and 
for the lOP are: adaptive phased-array radar, time-dependent system modeling, decon- 

volution and optimal control. 
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